############################## RCS限制性立方样条 #############################
##### 科研工具3  RCS限制性立方样条

library(readr)
library(corrplot)
library(rms)
setwd('父级路径')

demo <- read_csv("source.csv")


demoS2 <- demo
# demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')] <- log(demoS2[,c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')])

ddist<-datadist(demoS2)
options(datadist="ddist")
# RCSX <- c('MECPP','MnBP','MEHHP','MEOHP','MiBP','MiNP','MCOP','MCPP','MEP','MBzP')
tryCatch({
  alldataform
})
# RCSNK <- c(3,3,3,3,3,3,3,3,3,3)
tryCatch({
  RCSNKParam
})


i <- 1
Cairo::CairoTIFF(file="RCS.tiff", width=8, height=8,units="in",dpi=150)
# 默认25个 其他不能超过
par(mfrow=c(5,5))
for (i in 1:xsizeparam) {
  RCSX1 <- RCSX[c(-i)]
  paste0(RCSX1,collapse = '+')
  fml <- as.formula(paste0('Hyperuricemiaint ~ rcs(',RCSX[c(i)],',nk=',RCSNK[c(i)],')+',paste0(RCSX1,collapse = '+'),
                           '+Age + Sex + race + edu + FPL + activity + smoke + drink + Hypertension + CKD + urinecreatinine + BMI'))
  print(fml)
  fit <-lrm(fml,data=demoS2,x=FALSE)
  pred_OR<-Predict(fit,name=RCSX[c(i)],ref.zero=TRUE,fun=exp)
  
  ylim.bot<-min(pred_OR[,"lower"])
  ylim.top<-max(pred_OR[,"upper"])
  plot(pred_OR[,1],pred_OR[,"yhat"], 
       xlab = RCSX[c(i)],ylab = paste0("OR"),
       type = "l",ylim = c(ylim.bot,ylim.top),
       col="red",lwd=2)+
    lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)+
    lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)+
    lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) +
    points(1,1,pch=16,cex=1.2)#需要设置参考值,
}
dev.off()

